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Analysing borehole temperature data in terms of ground surface his- 

■ tory can add useful information to reconstructions of past climates. There- 
O ' fore, a rigorous assessment of uncertainties and error sources is a neces- 

' 5s ■ sar y prerequisite for the meaningful interpretation of such ground surface 

temperature histories. This study analyses the most prominent sources 
\ of uncertainty. The diffusive nature of the process makes the inversion 

Oh. relatively robust against incomplete knowledge of the thermal diffusivity. 

Similarly the influence of heat production is small. It turns out that for 
investigations of the last 1000 to 100 000 years the maximum depth of the 

■ temperature log is crucial. More than 3000 m are required for an optimal 
QQ ' inversion. Reconstructions of the last one or two millennia require only 
&\ , modestly deep logs (>300 m) but suffer severely from noisy data. 

m ■ 
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^ '■ 1 Introduction 

^P* , Borehole temperature data record past changes in ground surface temperature 

# ^ ' (GST) as the transient temperature signal diffuses into the subsurface. These 

^ , temperature data have been used to reconstruct time series of the ground surface 

$_i ' temperature in numerous studies. The value of such reconstructions as a proxy 

for the paleoclimate, h owever, is still subject of research and discussion (e.g . 



Mann & Schmidt 



2003 ). It has been shown (e.g. Gonzalez-Rouco et al. 20031 . 



Signorelli fc Kohll 12004 ) that surface air temperature (SAT) and soil tempera- 
ture correlate well at long timescales. Reconstructions have been carried out 
on various timescales from a few hundred to 100 000 years by several authors. 
While the shorter type serves, combined with proxy climate data, as a means 
to investigate natural and anthropogenic climate change, the longer reconstruc- 
tions study temperatures of the last ice-age and its transition to our current 
climate. 
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The problem is usually simplified to a purely conductive, one-dimensional 
one. If the model fails to capture reality, any inverse method will produce 
spurious results. Previous studies have already analyzed various effects that 
can possibly influence the transient signal. For instance, the error int roduced by 
the co mmo n assumptio n of one-dimensionality has been treated by IShen et al 



( 1995( l and iKohll ( 1999 ) . Advection of heat by groundwater flow can modify the 
signa l considerably ijClauser et al.lll997l . lkohil998llTaniguchi etlfll 19991 iRelter 
20051) . In the following analysis we will show that even with all simplifying 
assumptions fulfilled, it is possible to produce artefacts in the inverted time 
Although we use a particular algorithm for the study, the results will 



series. 



more or less apply to other programs as well. We will start with a short review 
of the theory of the algorithm in the next section. 



1.1 Theory 



The algorithms commonly used for GST history inversion assume a one-dimensional, 
purely conductive model. The physical proper t ies of the subsurface are known, 
the medi um can be layered ( Shen fe Beck!ll991 L 1992 ) or is assumed to be homo- 
geneous ( Mareschar&r^eltniml ^9^" Beltram^^^areschal 1995). Both ver- 
sions can be used on single boreholes or on ensembles of temperature logs in 
order to improve the signal to noise ratio. In a comparison of different codes 
I Beck et al.l Il992h . it was found that all of them perform equally well under 

similar circumstances. 

In the model considered here ( Mareschal fe Beltramilll992l . lBeltrami fe Mareschall 

19951) . the subsurface is considered homogeneous. The steady state temperature 
profile is defined by the heat flux at the surface qo, the pre-observational mean 
ground surface temperature To, thermal conductivity A, and heat production 
rate A of the subsurface. Together with a additiona l transient term T*(z,t\ 
the temperature at time t and depths z is given by IjCarslaw fc Jaegenll959l 
Beltrami fc Mareschallll995l ): 



T(z,t) = T Q + 



A 



Az 2 



+ T t (z,t). 



(1) 



In the particular case that the GST history can be parameterized by stepwise 



temperature changes T- 
becomes 



at times tj before present (t = 0) the transient term 



N 



Tt{z) 



3 = 1 



erfc 



crfc 



j'-i 



(2) 



If thermal conductivity, diffusivity and heat production are known, the only 
remaining unknowns are the steady-state GST, the surface heat flow and the 
values of the the GST history. If temperatures are recorded at depths z% i = 
1, M, a set of M linear equations can be formed from equation [H that can be 
solved for the unknown variables. The solution needs to be re gularized beca use 
the problem is ill-posed. The singular value decomposition ljLanczoslll96ll ) is 
used for this purpose by solving the linear problem and at the same time seeking 
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Figure 1: Time and depth scales of GST changes and respe ctive changes in to- 
day's subsurface temperature. Top: Idealized GST history l Haenel et al.|[l988l . 
table 10.9). The different line styles mark characteristic periods, respectively. 
Bottom: Temperature anomalies in the subsurface computed from the time se- 
ries. Line colours of the anomalies correspond to the periods of the upper panel. 
The dashed line is equivalent to the sum of the anomalies. 



a solution that minimizes the norm of the model. Implications of this procedure 
will be discussed in more detail in section [2J 

The original cod^j] was e mployed in various case studies for the recon- 
struction of past climates ( e.g. Beltrami et alJ 1997t Beltrami fc Bourlon 2004 , 
Clauser fc Mareschallll995l IClauser et al.lll997t Ijones et al.lll999h . It was also 
used in conjunction with other proxy data to improve the tempo ral resolution of 
the GST history as well as the coupling between GST and SAT {Beltrami et al 



1995llBeltrami fc Tavlorll995llSignorelli fc Kohj20ollNitoiu fc Beltramill2005h . 
Before analyzing the sensitivity of this algorithm it is worthwhile to investigate 
the effect of changing GST on the temperature measured in boreholes. 



1.2 Paleoclimate influence on subsurface temperatures 

The magnitude and time/depth scale of GST changes can be seen in figure [TJ 
We used a temper ature history of the last 1 million years based proxy data 
l Haenel et al. 19881 . table 10.9) to demonstrate the effect of this transient signal 
on the subsurface temperatures. The curve is valid only for Switzerland but it 
displays the main features of the climate of the last 1 million years. The proxy 
curve is subdivided into 3 parts and for each a perturbation with respect to the 
mean is calculated. The transient signals for each of the periods are computed 



1 A MATLAB version of the program including a GUI and the optimal choice of 
the regularization parameter by means of the L-curve method can be downloaded from: 
http : //uuu . geophysik . ruth-aachen . de/Dounloads/sof tuare . html 



separately (figure[IJ lower panel, solid lines). The mean GST (To in equation[TJ 
for each of the computations is assumed to be 0°C. 

The total transient signal (dashed line) is then the sum of the three partial 
signals. It is apparent that changes of the order of 10 6 years for this particular 
GST history will appear in a temperature log as a virtually straight line with 
a slope of about -1 K km" 1 unless the temperature log runs down to a depth 
of at least 4000 m. The slope will change if a different GST history is used 
but it is clear that this part of the transient signal cannot be resolved by an 
inversion, especially not in the presence of noise. Thus, fluctuations of the pe- 
riod 100 000 to 1000 000 years will usually appear only in the pre-observational 
mean (POM). The glacial temperature minimum and the postglacial temper- 
ature increase (100 000 to 10 000 years) cause the most prominent disturbance 
of the temperature profile. The temperature deviation reaches its maximum 
amplitude at about 1500 m depth, but still has a value of 1.5 K at 4000 m 
depth. Because most of the available undisturbed temperature data are from 
smaller depths, a large portion of this signal is frequently not recorded. The 
signal due to temperature variations during the Holocene can be traced down to 
a maximum depth of 1000 m. Although temperature data in this depth range 
are more numerous, reconstructions have been attempted using substantially 
shallower data. 

From the discussion above two questions arise. Section [2] focuses on the 
question what the value of the regularizing parameter e for any given problem 
should be and whether some objective criterion can be used to define it. In sec- 
tion [3] the question is addressed if a reliable reconstruction of past temperatures 
is possible when only part of the transient signal is contained in the data. 



2 Optimum choice of regularization parameter 

It was mentioned in the introduction that a regularization is needed to solve the 
inverse problem in a stable m anner. For th e singular value decomposition, two 
methods are commonly used ( Menke 1989h . Either only the p largest singular 



values are used in the inversion or a damping parameter e is added to the singular 
values such that 

1 Aj 

The stabilization of the solution is obtained at the cost of loss of model resolu- 
tion. Both methods require a choice of a regularization parameter /?, i.e. the 
number of singular values discarded, or the magnitude of the damping parame- 
ter. A regularization of this type tends to minimize the norm of the model. This 



is equivalent to minimizing the objec tive function </>(m) given by (jAster et al 



is equivalent to minimizing the objec ti 
2004L iFarquharson fc Qldenburgll200l : 



0(m) = ||G(m)-d||+^||m||. (4) 

Here the first term represents the L2-norm of the data misfit, with G(m) being 
the forward model and d being the data. The second term represents the L2- 
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increasing e 



Model norm ||m|| 



Figure 2: Principle of the L-curve. With increasing e the model norm decreases 
and the data norm increases. The optimum value for e can be found at the 
maximum curvature point of the curve. 



norm of the model m. In terms of paleoclimate inversion, minimizing the model 
norm can possibly change the amplitude of paleoclimatic disturbances. Too 
small an e, on the other hand, will lead to spurious oscillations in the solution 
of the GST. Thus it is important to choose a proper value for the damping 
parameter, preferably based on some objective criterion. 

One way to gain insight in to the problem is the L-curve ( Hansen 19981 . 
Farquharson fc Oldenburg! 12004 ). so called for its shape. To create this curve, 
the inversion is repeated for a range of regularization parameters and the data 
norm is plotted versus the model norm (figure || on a logarithmic scale. For 
large values of e the model norm is small but the data norm is large. A small 
e leads to a large model norm and to a data misfit norm that is close to the 
variance of the data. These two cases define the asymptotic behaviour of the 
L-curve. At the maximum curvature point o f the curve, a compromise is found 
between the magnitude of those two norms ( Hansen 1998t ). This is the optimal 
value for the regularization parameter. 

We studied our inversion problem with the help of the L-curve criterion. A 
synthetic temperature log was created using the parameter set given in table [TJ 
The transient temperature anomaly given in figure Q] was added to the steady 
state solution and random, normally distributed noise was added. The noise 
amplitudes for the test data sets were 0.01, 0.2, 0.4, and 0.6 K, respectively. 

Figure [3] shows the variation of the L-curve and the optimum e as a function 
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Table 1: Parameters of the synthetic model used for the L-curve studies. 



Parameter 


Value 


Unit 


Sampling rate 


20 


m 


Total depth of log 


5000 


m 


Number of time steps for inversion 


20 




Start time of inversion 


1000 


ky b.p. 


Stop time of inversion 


0.1 


ky b.p. 
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Figure 3: L-Curve for different noise levels. 
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Figure 4: Inverted GST-history for different noise levels and optimal e from 
figure [3j The curves are each shifted by 5°C for better visibility. 

of the noise level. It is apparent that only for the smallest noise level the L- 
curve takes on a well shaped form with an maximum curvature point easy to 
determine. For all noise levels the L-curve displays a straight line region where 
any change in e leads to large changes in the model norm but does not change 
the data fit at all. This means that a value of e less than the value found in the 
maximum curvature point is not justified by the data. It is noteworthy that the 
model norm of the optimal model does not vary much with the increasing noise 
level. This indicates that the inverted model is robust against noise for this 
configuration. This point is illustrated by comparing the inverted GST histories 
for optimal values of e (figure [4]). The major feature is the temperature increase 
at the end of the last ice age. It is preserved in all of the results in comparable 
magnitude and robust against noise. Smaller temperature variations during the 
Holocene period, on the other hand, are only resolved at the smallest noise level. 
The L-curve provides a tool to answer the important question if a sought for 
feature in the paleoclimate record can be resolved by the data. 

Finally, figure [5] shows a comparison of inversion results for a given noise 
level and a range of e-values. The optimal value would be about 0.14, yielding a 
rather smooth model. The interpreter, possibly aware of a paleoclimatic proxy 
curve like the ones shown in the figure, might be tempted to choose a smaller 
value that seems to better fit the temperature variations in the model during 
the Holocene. However, figure [3J indicates that the data fit does not improve at 
all by choosing an e smaller than 0.14. The small amplitude variations on top of 
the temperature increase at the end of the ice age appear to be only numerical 
artifacts. 
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Figure 5: GST-histories for a fixed noise level of 0.2 K and varying values of 
the damping parameter e. The optimal choice of e according to the L-curve 
criterion is 0.14. 

3 Sensitivity analysis 

3.1 Length of the temperature log 

In the introduction the question was posed if the full depth range of the transient 
signal is needed for a valid reconstruction of the GST history. To analyze this 
problem, a synthetic temperature log with a simplified GST history of the last 
1 million years is computed. The parameters for the temperature log are taken 
from table [21 The GST history that produces the transient temperature signal 
(grey line, figure [6]) does not represent an actual paleoclimatic curve but rather 
represents an idealized history with correct order of magnitude for amplitude 
and timing of the events. 

For the inversion procedure, the originally 3000 m deep log was cut off at 
depths of 2000 m and 1500 m. All three logs were inverted with and without 
synthetic noise applied. We chose normally distributed noise with a standard 
deviation of 0.3 K. For the inversion of the noise free and noisy logs, a constant 
damping parameter e of 0.1 and 0.2, respectively, was used. 

The resulting GST histories for the different inversion runs are displayed in 
figure [6l For both the noisy and the noise free data a log depth of 3000 m seems 
to be sufficient to fully recover the amplitude of the last glacial stage. A decrease 
of the maximum log depth also reduces the amplitude of the glacial cooling. 
Further the phase of the events is shifted to later times. The effect is more 
pronounced for the noisy data, but even for noise free data any reconstruction 
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Figure 6: A shortening of logs results in different GST reconstructions. Solid 
line represents the GSTH that was used to drive the forward model. The dashed 
lines show reconstructions for different depths of the logs. The temperature of 
the last glacial minimum is only resolved for the 3000 m log. Not even the 
timing of the climatic optimum in the Holocene is correctly resolved for the 
1500 m log. 



Table 2: Values for the thermal properties that are used throughout the text to 
compute synthetic temperature-depth profiles. 



Parameter 


Value 


Unit 


Thermal conductivity 


2.5 


W (m K)" 1 


Volumetric heat capacity 


2.30 -10 6 


J m 3 K 1 


Thermal diffusivity 


1.09 -10- 6 


2 -1 

m s 


Heat production rate 


1.0 -10- 6 


W m" 3 


Surface temperature 


10.0 


°C 


Heat flow 


0.06 


W m 2 


Sampling rate 


1.0 


m 
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Figure 7: Upper panel: Various reconstructions of paleo temperatures of the 
last two millennia. Lower panel: Inversio n results for a 2000 m deep noise free 
temperature log using the proxy-curve of Esper et al. ( 20021 ) as input data. 



of the last ice age paleoclimate based on logs shallower than 2000 m seems to 
be problematic. 

A similar procedure of shortening logs can be used to estimate the necessary 
depth of temperature logs for millennial scale climate reconstructions. Cur- 
rently plenty of proxy reconstructions for the last one or two millennia exist and 
their validity is debated. The geothermal method can provide valuable insight 
into this problem as it is the only direct measurement of paleo-temperatures. 
The more it is important to analyze the information contained in the tempera- 
ture data. Figure [3 upper panel, shows a few reconstructions of the Northern 
hemisphere annual mean SAT together with a reconstruct i on based on boreho le 



temperature data ijMann fc Jonesll2003l . lEsper et al.ll200l iHuang et alllioooh 



mg 

In our analysis we use the reconstruction of lEsper et alJ ([2002]) as the forcing 



to compute the transient temperature perturbation. The other parameters of the 
forward model are taken again from table [21 Additionally, a pre-observational 
mean needs to be specified. As the proxy curve contains no information about 
this parameter, it was arbitrarily chosen to be the mean temperature in the 
interval 800-1200 years before present. FigureO lower panel, displays an inver- 
sion result for this synthetic temperature log. The temperature data are noise 
free and the log runs 2000 m deep. Because of the favourable conditions for 
this particular inversion run the information contained in this reconstruction 



10 



can be considered the maximum knowledge that can be extracted by inverting 
temperature-depth data. The warming trend of the last 200 years is captured 
very well. No insight into climate variations of the preceding centuries can be 
gained, but the pre-observational mean is preserved in the reconstruction. Thus 
the borehole method effectively provides two unknown pieces of information: 
The temperature change from the average temperature of the last millennium 
and the mean temperature before that period. 
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Figure 8: Inversion results for decreasing maximum depth of the temperature 
log. Temperature logs were cut off at various depths. Thermal parameters are 
the same as in table El Upper panel: Input temperature data for the inversion 
was noise-free. The damping parameter e equals 0.1. Lower panel: Temperature 
data are disturbed by normally distributed noise with a standard deviation of 
a = 0.1 K. A higher e of 0.6 is used in the inversion. 

These considerations are only valid for the best available data. In reality, 
data are noisy and boreholes are not sufficiently deep to allow extraction of this 
information from the temperature data. Figure [8] compares inversion results 
where only shorter logs are available and noise is added to the data. For the 
case without noise, the average temperature of the last millennium can be re- 
solved for logs that run deeper than 300 to 500 m. For shallower logs a strong 
decrease in amplitude of the GST reconstruction is apparent. If noise is present 
in the data even the 1000 m temperature is not able to fully recover the mean 
temperature of the last 1000 years. These observations are different from the 
results for glacial/post-glacial reconstructions. The reason is that the amplitude 
of the transient temperature signal due to temperature changes in the last few 
thousand years is one order of magnitude less than that observed for the post-ice 
age temperature rise (see figure [J). The input signal used here is a global SAT 
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Figure 9: Estimating the transient temperature anomaly and undisturbed gra- 
dient results in a biased estimate for the gradient and a smaller temperature 
anomaly. 



reconstruction. At a specific site the amplitude of the transient signal can be 
higher and thus improve the signal to noise ratio that we assume in our example. 
Thus, in practice the expected signal should be compared to the noise present 
in the data to assess if the desired signal can be extracted from the data. 

Both examples show a strong dependence of the inverted GST history on 
data quality and both fail to reproduce the input forcing under certain circum- 
stances. The ill-posed nature of the problem together with the regularization 
using the damping parameter e cause this behaviour. Besides preventing numer- 
ical instability in the inversion the procedure also minimizes the norm of the 
parameter vec tor. Usually this effe ct is desired to avoid unnecessary compli- 



cated models ([Constable et al.lll987P . However, in our case a minimal solution 



in this sense will have smaller GST variations compared to the GST variations 
that created the transient temperature anomaly. This is especially true since 
the mean GST To and the surface heat flow go are also solved for. Figure [9] illus- 
trates this problem. The black line denotes the true temperature perturbation 
due to a step increase in ground surface temperature. If one inverts simultane- 
ously for GST history, qo, and To using a minimum norm model in the presence 
of noise, a steady state temperature profile (grey line) tilted toward the temper- 
ature anomaly will be obtained instead the true one (dashed black line). In a 
case with a step increase in GST, the estimate for q will be too low and that for 
To will be too high. An example of this trend will be given in the case example. 
It is also notable that the temperature perturbation crosses the steady-state 
profile at two points (arrows). Thus, a false GST history is created with a de- 
crease in ground surface temperature and a subsequent increase although the 
input model uses only a step increase. Even with an e chosen optimal for the 
inversion, this problem will persist as it is imminent in the regularization. 
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Figure 10: Systematic changes in thermal diffusivity cause shifts in the timing 
of the events. Thermal conductivity is varied from 70 % to 130 % of the original 
value. For better visibility, inverted curves are stacked and plotted as smooth 
lines rather than a series of step functions. 



3.2 Thermal parameters 

The inversion procedure requires that thermal conductivity, thermal diffusivity, 
and heat production rate are exactly known. In a real situation errors for this 
parameters can be minimized but cannot be avoided altogether. For instance 
errors might be introduced by preferential sampling of units that are not rep- 
resentative of the whole sequence. Also, missing or false corrections for in-situ 
conditions of temperature and pressure can cause systematic deviations from the 
correct thermal parameters. In the following the uncertainties in the inversion 
resulting from these errors are considered. 

In general, the inversion procedure interprets curvature in the temperature 
log as a transient signal. Thus, only the terms introducing curvature in equa- 
tion Q] will have an influence on the inversion result. The transient term of the 
model has thermal diffusivity as its only relevant parameter. We can expect 
that this property will have the largest influence on the solution. Within the 
steady state terms, only the heat production term introduces curvature. As 
heat production rate is usually small, the influence of this term is expected to 
be less important for boreholes of moderate depth (see I3-2|) . Heat production 
and thermal conductivity are related inverse proportional in this term. Thus, 
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Figure 11: Systematic changes in heat production rate cause variations in the 
amplitude of the events. Heat production rate is varied from 70 % to 130 % of 
the original value. The inverted curves are plotted as smooth lines rather than 
a series of step functions for better visibility of the effect. 
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to analyze the errors introduced by this term, only on of the parameters needs 
to be studied. We chose the heat production rate. Thermal conductivity has its 
main influence in the steady-state terms for heat conduction. A false assumption 
will therefore lead directly to a proportional error in the heat flow density but 
does not affect the transient part of the temperature log. In practice, of course, 
thermal diffusivity is indirectly determined as the ratio of thermal conductivity 
and volumetric heat capacity and will have an indirect effect on the solution. 

In order to support this qualitative discussion with quantitative results a 
synthetic temperature log was inverted with false assumptions of thermal diffu- 
sivity and heat production rate in two separate inversion runs. The parameters 
were varied by ±30% of the original value. A simplified GST history was used 
in the forward model to compute the transient temperature anomaly. The GST 
history used is comprised of several step changes that aim to reproduce the 
general magnitude and timing of the expected effects and do not represent an 
actual GST history. The results of these experiments are shown in figures [10] 
and [TTJ Thermal diffusivity influences the timing of events in the reconstructed 
time series. A phase shift of the entire reconstruction occurs when a false es- 
timate is used in the inversion. The amplitude of computed past temperature 
remains unchanged. Heat production influences the magnitude of the paleocli- 
matic events but retains timing of events. The effect of the heat production 
increases with time, consistent with the fact that the heat production term in- 
creases quadratically with depth. It is apparent, that the heat production will 
only be significant for very large deviations from the correct value and deep 
boreholes. For instance, it should be irrelevant for millennial scale studies. 
Thermal diffusivity can change the timing of events on all timescales and it is 
also important for inversions of much shorter timescales than considered here. 



3.3 Heterogeneity 

The previous discussion considered a homogeneous subsurface only. The ques- 
tion is of course whether the conclusions drawn a valid for a heterogeneous 
medium. To address this question we will assess the impact of a heterogeneous 
ID layered earth model on the results of the in version. In the followin g the 



transient profiles computed by a FD algorithm IjRath fc Mottaghyl 120051 ) will 
be considered as ground truth. Random variations of the thermal properties 
on distinct length scales will be introduced and the inversion results studied. 
We will discuss the effect of systematic variations on a synthetic example of a 
sedimentary rock column. We will study the steady state T-z profile that results 
from compaction and elevated temperatures. 

For the forward computation a simplified GST history of the last 100 000 
years is used, with a pre-observational mean of 5°C. The profile is computed 
down to 10 000 m with a vertical resolution of 5m. For the analysis only the top 
3000 m are used, the deeper parts are only necessary to avoid edge effects in 
the solution. Thermal parameters are those of table [2] if not explicitly specified. 

Obviously, the largest influence on the temperature profile will come from 
variations in thermal conductivity. Deviations from the straight line that is 
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assumed for the stationary temperature will be falsely interpreted as noise or 
paleoclimate. In the homogeneous subsurface framework it is impossible to 
study variations of A and K separately. One way to overcome this problem is to 
cast the steady state part of the equation in terms of the Bullard depth Zb, for 
the homogeneous case defined by z& = z/X. Consider the differential equation 
for ID heat transport in a homogeneous medium: 

d 2 T _ldT A 

Substituting dz = Xdzb yields 

d 2 T 1 dT A 

Thus if k' = ft/ A 2 and A' = A/X 2 , equation Q] can be written as: 

A'z 2 

T(z b ) = T Q + q z b H — ^ 



Using the transformed variables and taking the heat production term to the left 
hand side, the steady state part of the profile can be reduced to a homogeneous 
problem. With this construct it is possible to use a layered model for the thermal 
conductivity and a homogeneous one for the thermal diffusivity. This allows 
us to study the influence of thermal diffusivity alone in the algorithm. With 
the general framework established, we will now turn to the different aspects of 
heterogeneity. 

The noise on varying length scales has been studied by modifying the FD 
model to compute the temperature profiles. The layers of the model were as- 
signed random values of thermal conductivity with a mean of 2.5 W(m K)" 1 
and a standard deviation of 0.5 W(m K)" 1 . Volumetric heat capacity was held 
constant throughout the experiment, resulting in a simultaneous variation of 
thermal conductivity and diffusivity. Ten different configurations were consid- 
ered with layer widths x increasing from 5 m to 1000 m. For each configuration 
50 realizations of the thermal conductivity-depth profile were generated. These 
were then used in the forward model to compute 50 temperature depth profiles. 
Additionally, temperature profile for an equivalent half space and for a layered 
representation (equation [8]) were computed. The differences between the "true" 
temperature profile and the two profiles with assumptions are shown in figurefT2l 
The misfit of the equivalent homogeneous model is large because variations in 
the steady state temperature profile are not reflected in the homogeneous model. 
This highlights the importance of the knowledge of thermal conductivity. On 
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Figure 12: Misfit between temperature profiles computed from the FD-model 
and homogeneous models for varying length scales of variability. Mean RMS- 
values are given for a homogeneous subsurface (o) and for homogeneous k only 
(o). Grey areas denote min/max range. 



the other hand, for the model with homogeneous k only deviations of the order 
of a few mK occur, suggesting that random variations of thermal diffusivity are 
smoothed out. 

Next, the random data have been used as input to the inversion procedure. 
For each inversion the value and timing of the maximum temperature in the 
Holocene have been recorded in order to compress the results (figure [13]) . To 
accommodate the increasing misfit between model and data for increasing noise 
length scale, the damping parameter e was increased from 0.35 to 1 according 
to an increase in noise length scale from 5 m to 1000 m. This leads to a stronger 
damping and is reflected in the decrease of the maximum temperature value 
because models with smaller norms are preferred for large e. The influence of 
noise on the result becomes most pronounced for length scales larger than 100 
m when the scale of the signal is on the order of the noise scale. The timing of 
the event is modified in a similar way. 

In crystalline rocks random variations of thermal properties might prevail, 
in sedimentary rocks systematic changes are much more important. Changing 
depositional environments throughout time create diverse rock types. The in- 
fluence on the inversion depends on the particular setting. One aspect that is 
often neglected is the variation of thermal properties in uniform sequences of 
rock due to the temperature dependence of thermal properties and the com- 
paction with increased overburden. In figure [14] a geologic profile of a single 
rock type is considered. The thermal conductivity of the matrix is taken as 3.5 
W(m K)" 1 . An exponential decrease of porosity due to compaction is assumed 
(A thv's law) with a su rface porosity of 0.45 and a characteristic depth of 750 
m ((Alien fc Allenlll990h . Thermal conductivity of the matrix is taken inversely 
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Figure 13: Variation of the timing and amplitude of the maximum temperature 
in the Holocene versus the length scale of the random noise as defined in the 
text. Left: Value of the maximum temperature. Right: Timing of the maximum 
value. Grey area denotes 75 % quartiles of all realizations. 
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Figure 14: Systematic variation of thermal properties with depth, a) Poros- 
ity, temperature, and reduced temperature; b) thermal conductivity and heat 
capacity of the rock matrix, depending on temperature; c) effective thermal 
conductivity and diffusivity with porosities from panel a); d) Effective thermal 
properties for constant porosity of 0.1. 
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proportional to temperature ( Haenel et al. 1988f l and heat capacity of the ma- 
trix is approximately proportional to temperature in the range considered. For 
simplicity the pore fluid is assumed to be fresh water, thermal proper ties are 
taken from reference tables l|Wagner fc Pruf^OOl iKaye fc Laby)ll968f l. Tem- 
perature and heat flow at the surface are 10 °C and 0.06 W uf 2 , respectively. 
These assumptions lead to a curvature in the thermal properties (figure [14] c) . 
The curvature is reflected in the reduced temperature T r (figure [14] a) . This 
effect is similar in magnitude to the influence of paleoclimate. It is apparent 
that these effects need to be taken into account for any paleoclimate inversion. 

For comparison, effective thermal rock properties are also given for a constant 
porosity of 0.1 to consider only the temperature/pressure effects (figure fl4l d). 
The shape of the A — z profile changes considerably displaying the prominent 
influence of the porosity decrease. The k — z profiles show that the porosity 
effect is compensating the temperature effect. Diffusivity varies by about 20 % 
in this example. When comparing these results to the inversions of the previous 
sections it is clear that thermal conductivity needs to be known accurately to 
achieve reliable results. Diffusivity is again somewhat less important. 



4 Case study: KTB 



A temperature log assumed to be in equilibrium was record ed in the KTB pilot 
hole down to 3990 m on September 17, 1997 and analyzed in Clauser ( 19991 ) . An 
earlier recording, dating from February 1996, was analyz ed in detail for ther - 



mal processes, including transient effects of paleoclimate l|Clauser et alJll997n 



Shallow temperature logs in the vicinity of the KTB site have also been inter- 
prete d in terms of paleoclimatic influence ( Clauser fc Mareschall 1995 . Clauser 
1999l l. In the light of the previous sections, we will revisit the analysis of the 
KTB borehole and discuss the uncertainties connected with the interpretation. 
In that study the temperature data was inverted for paleoclimate in the period 
from 10 5 to 10 2 years before present, using 100 time steps, a thermal conductiv- 
ity of 2.92 W (m K)" 1 , a thermal diffusivity of 10~ 6 m 2 s" 1 , and a heat production 
rate of 1.1 /iW m 3 . Their analysis yielded a peak to peak amplitude of 10 K 
for the temperature increase from the latest glacial stage to the Holocene. 

The geological profile in the KTB pilot hole is primarily composed from two 
lithologies: Metabasites and Gneisses. These can easily be separated using the 
Gamma-ray- log (SGR) (figure [HI left panel). Thermal conductivity (figure HH 
right panel) does not seem to vary systematically over the depth of the bore- 
hole, justifying the assumption of a homogeneous half-space. The mean thermal 
conductivity is (2.93 ± 0.60) W(m K)" 1 . This value reduce s to 2.82 W(m K 



if a temperature-pressure correction is applied to the data (|Buntebarthlll99l[ l 



V( 
th 



Knowledge of the thermal diffusivity is much more uncertain. Literature cites 
ranges for Gneisses and Am phibolites as 0.5 to 1.2 - 10~ 6 m 2 s" 1 an d 0.6 to 
0.8 • 10~ 6 m 2 s 1 , respectively |Landolt-Bornsteinl[l98l . ISeipoldl <|l995h reports 
(0.8 ±0.2) • 10 -6 m 2 s" 1 for Amphibolite samples from neighbouri ng sites. Again , 
a correction for temperature and pressure of the value given by lSeipoldl (<1995l ) 
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Figure 15: Composite log, KTB pilot hole. Panels from left to right: 1) Spec- 
tral Gamma Ray (SGR), temperature (TEMP), temperature gradient (TMPG); 
2) Generalized lithology; 3) Thermal conductivity: Core measurements (TC), 
mean value (mean TC), ±lcr of the mean value (shaded area), corrected for tem- 
perature and pressure (corr TC); 4) Heat production A from spectral gamma 
ray (A), mean value of A (mean A), thermal diffusivity re (TD). 
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Figure 16: L-curve for the inversion of the KTB pilot hole-temperature log. 

reduces the mean diffusivity to a value of 0.70 m 2 s" 1 . 

The next step in the analysis is the choice of the optimum regularization 
parameter e. Figure fl6l shows the L-curve for varying values of e. The optimum 
choice seems to be e = 0.4. FigurefTTI highlights this in more detail. The inverted 
time series for the different values of e are shown. The direct comparison of 
results is especially helpful to analyze which features of the GST history are 
numerical oscillations and which ones represent actual information. At a level 
of 0.4 oscillations that change the warming amplitude seem to be dampened 
out. A more conservative interpreter might choose an even higher value of 0.6, 
but in this range of e-values the main features are only changed by a few tenth 
of a degree Celsius. 

Following the choice of the damping parameter, we shortened the log to illus- 
trate the detrimental effect of the log length. The maximum depth of originally 
3990 m was successively cut to lengths of 3000, 2000, and 1500 m. Figure fTHl 
shows the reconstructed GSTH and transient temperature perturbation in the 
borehole for various cut off depths. Deviations from the expected results be- 
come significant for depths less than 2500 m. It is noteworthy that not only 
the temperature drop of the last glacial stage decreases but also the timing of 
the optimum temperature in the Holocene changes. An interpreter of paleocli- 
matic inversions must be aware that the post glacial temperature increase will 
influence the shape and timing of later events as well. 

Table [3] shows the inversion results for the undisturbed GST (i. e. the pre- 
observational mean), the heat flow into the model and the undisturbed gradient. 
This highlights that the systematic variation of the GST history is accompanied 
by an according change in the steady-state parameters. As discussed in the 
section on length of temperature logs, the heat flow value decreases whereas 
the GST increases. This is consistent with inverting a step increase in ground 
surface temperature using a log that is too short too fully resolve the GST 
history. It can be expected from this data that crustal heat flow estimates may 
be too low by about 8 mW nf 2 , if a transient signal of the last ice age has not 
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Figure 17: Inverted GST histories for varying values of the damping parameter. 




Figure 18: Inversion results for the KTB log for different depths of the log. 
Top: Inverted GST history for different log length. For clarity, the step ampli- 
tudes of the GST histories are shown as continuous lines. Bottom: Transient 
temperature signal in the subsurface. 
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Table 3: Inversion results for the KTB log, cut off at various depths. The 
inverted parameters vary systematically with the depth of the log. 



Depth (m) 


To (°C) 


q (mW nT 2 ) 


dT/dz (mK nf 1 ) 


4000 


4.14 


86.3 


29.6 


3000 


4.70 


84.8 


29.6 


2000 


5.51 


81.8 


28.0 


1500 


6.15 


78.1 


26.7 



been properly removed. 

5 Conclusions 

The sensitivity analyzes in the previous sections have shown that the inversion 
results are strongly dependent on data quality. The parameters most important 
depend on the time scale of inversion. For inferences about the post-glacial 
temperature rise the maximum depth of the temperature log is the most chal- 
lenging requirement. If it is met, the inversion is robust against noise in the 
temperature data. A depth of more than 3000 m seems to be optimal to re- 
solve the temperature rise adequately. Unfortunately, undisturbed temperature 
logs running that deep are rare. Reconstructions of the last millennium require 
only modestly deep logs. A depth of about 300 m to 500 m is suitable to de- 
rive a mean temperature for the last 1000 years. The exact values will depend 
on the magnitude of the noise and the transient signal at a specific site. For 
instance the IHFC database for borehole temperatures and climate reconstruc- 
tion lists 473 out of 754 temperature logs with a maximum depth greater than 
300 m (http: //www. geo.lsa.umich.edu/~climate). However, for this type of 
inversion the demand on data fidelity is very high. Even a noise level of 0.1 K 
seriously degrades the inversion results. In practice the noise level a 

Errors in the thermal properties, namely thermal diffusivity and heat pro- 
duction, seem not to be the main source of error. The diffusive nature of the 
signal and the spreading over a wide depth range greatly dampen the influ- 
ence of the thermal diffusivity. The importance of heat production increases 
with reconstruction length. But even for studies of the last 100 000 years er- 
rors introduced by false assumptions of the heat production rate will be small 
compared to other sources of error. 

The studies on the regularization parameter show that a thorough noise 
analysis should accompany any attempt to invert GST histories from geother- 
mal data. The correct choice of the regularization is crucial and should be 
justified by some objective means. A non-optimal choice of the regularization 
can easily lead to artefacts that are much larger than the calculated model 
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errors. Especially to small an e introduces events that could be mistaken for 
variations of paleo-temperatures. We used t he L-curye crit e rion but other meth 
ods l i ke generalized cro ss validation (e. g., Hansenl 19981 lHaber fc Oldenbur 



Hansen! Il99 



2000L lAster et al. 20041) or truncated iterative methods (e 
Hankdll995l lBerglundll2o"o2h could also be utilized. 

Finally, although the least-squares fit with a Tikhonov regularization mini- 
mizing the model norm is a generally applied method in solving ill-pos ed prob- 
lems, it might not necessarily be the best choice, as pointed out by iHansenl 
1 19981 ). Other regularization methods might better preserve the information 
contained in borehole temperatures . These include Bayesian-type smoothing ap- 
proaches JSerban fc Jacobsenl 2001 ) as well as nonlinear minimum support reg- 
ularizers (jPortniaguine fc Zhdanov 19991 . Zhdanov 20021 . Zhdanov fc Tolstaval 
2004T I . Research in the development of novel approaches to regularization and 



the incorporation of prior knowledge are an urgent task for the future. 
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